{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.stats import multivariate_normal\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "import bayesnet as bn\n",
    "\n",
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARUAAAD8CAYAAABZ0jAcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD2RJREFUeJzt3X+sZGV9x/H3p+uSbLekW+X36qaakLVrVMAbpJZUqApC\nagHTJpAGjTHZ2KipptkEYmL9z9aNNbFR7NqSalIxbWSR6OKWNW2oNVh3+bG7CKuoELms/NKFqhsF\n+u0fcy4Ol/tj7s5zZ+7cfb+SyT3nOc9z73Nm2A9z5pwz31QVktTKb4x7ApJWF0NFUlOGiqSmDBVJ\nTRkqkpoyVCQ11SRUklyf5NEkB+fZniSfTHJ/kv1Jzunb9tYkh7pt17SYj6TxafVO5Z+Bty6w/RLg\nzO6xFbgOIMka4FPd9i3AVUm2NJqTpDFoEipVdRvwkwW6XAZ8vnpuBzYkOR04F7i/qn5QVb8Cvtj1\nlTShXjSiv7MR+FHf+kNd21ztr5/rFyTZSu9dDuvXr3/dK1/5yuWZqST27dv3eFWdfCxjRxUqQ6uq\nHcAOgKmpqdq7d++YZyStXkkePNaxowqVaeBlfesv7drWztMuaUKN6pTyzcA7urNA5wFPVtVh4NvA\nmUlenuQE4Mqur6QJ1eSdSpIbgAuAk5I8BPw1vXchVNVngF3ApcD9wC+Ad3XbnknyPmA3sAa4vqru\naTEnSePRJFSq6qpFthfw3nm27aIXOpJWAa+oldSUoSKpKUNFUlOGiqSmDBVJTRkqkpoyVCQ1ZahI\naspQkdSUoSKpKUNFUlOGiqSmDBVJTRkqkpoyVCQ1ZahIaspQkdSUoSKpqVZlTxcsXZpkW5K7usfB\nJM8meXG37YEkB7pt1t2QJtzQ31HbV7r0LfSKgX07yc1V9Z2ZPlW1Hdje9X8b8MGq6q9oeGFVPT7s\nXCSNX4t3KkstXXoVcEODvytpBWoRKvOVNH2BJL9Jr5D7l/qaC9iTZF9X2lTSBBt12dO3Af8969Dn\n/KqaTnIKcGuS+7qC78/TX0t506ZNo5mtpCVr8U5lvpKmc7mSWYc+VTXd/XwU2EnvcOoFqmpHVU1V\n1dTJJx9T3WhJI9AiVAYqXZrkt4E3Al/ua1uf5MSZZeAi4GCDOUkak6EPf+YrXZrkPd32z3RdrwD+\nvap+3jf8VGBnkpm5fKGqvjbsnCSNT3oVSSfL1NRU7d3rJS3Sckmyr6qmjmWsV9RKaspQkdSUoSKp\nKUNFUlOGiqSmDBVJTRkqkpoyVCQ1ZahIaspQkdSUoSKpKUNFUlOGiqSmDBVJTRkqkpoyVCQ1ZahI\naspQkdSUoSKpqVHVUr4gyZN99ZQ/POhYSZNlJLWUO/9VVX98jGMlTYhx1FJuNVbSCjTKWspvSLI/\nyS1JXrXEsSTZmmRvkr2PPfZYg2lLWg6j+qD2DmBTVb0G+HvgpqX+AsueSpNhJLWUq+qpqvpZt7wL\nWJvkpEHGSposI6mlnOS0dLVNk5zb/d0nBhkrabKMqpbynwJ/keQZ4ChwZfXqrc45dtg5SRofaylL\negFrKUtaMQwVSU0ZKpKaMlQkNWWoSGrKUJHUlKEiqSlDRVJThoqkpgwVSU0ZKpKaMlQkNWWoSGrK\nUJHUlKEiqSlDRVJThoqkpgwVSU2Nquzpn3c1fw4k+WaS1/Zte6BrvyuJ3xEpTbhRlT39IfDGqvpp\nkkuAHcDr+7ZfWFWPDzsXSeM3krKnVfXNqvppt3o7vfo+klahUZY9nfFu4Ja+9QL2JNmXZOt8gyx7\nKk2GoQ9/liLJhfRC5fy+5vOrajrJKcCtSe6rqttmj62qHfQOm5iampq8uiLScWIkZU8BkrwG+Efg\nsqp6Yqa9qqa7n48CO+kdTkmaUKMqe7oJuBG4uqq+29e+PsmJM8vARcDBBnOSNCajKnv6YeAlwKe7\nksrPdNXPTgV2dm0vAr5QVV8bdk6Sxseyp5JewLKnklYMQ0VSU4aKpKYMFUlNGSqSmjJUJDVlqEhq\nylCR1JShIqkpQ0VSU4aKpKYMFUlNGSqSmjJUJDVlqEhqylCR1JShIqkpQ0VSU6Mqe5okn+y2709y\nzqBjJU2WoUOlr+zpJcAW4KokW2Z1uwQ4s3tsBa5bwlhJE2QkZU+79c9Xz+3AhiSnDzhW0gQZVdnT\n+foMXDLVsqfSZJiYD2qrakdVTVXV1Mknnzzu6UiaR4tayoOUPZ2vz9oBxkqaICMpe9qtv6M7C3Qe\n8GRVHR5wrKQJMqqyp7uAS4H7gV8A71po7LBzkjQ+lj2V9AKWPZW0YhgqkpoyVCQ11eKUsjR2N905\nzfbdh3j4yFHO2LCObRdv5vKz57yOUsvMUNHEu+nOaa698QBHn34WgOkjR7n2xgMABssYePijibd9\n96HnAmXG0aefZfvuQ2Oa0fHNUNHEe/jI0SW1a3kZKpp4Z2xYt6R2LS9DRRNv28WbWbd2zfPa1q1d\nw7aLN49pRsc3P6jVxJv5MNazPyuDoaJV4fKzNxoiK4SHP5KaMlQkNWWoSGrKUJHUlKEiqSlDRVJT\nhoqkpoYKlSQvTnJrku91P39njj4vS/IfSb6T5J4kf9m37SNJppPc1T0uHWY+ksZv2Hcq1wBfr6oz\nga9367M9A/xVVW0BzgPeO6u06Seq6qzusWvI+Ugas2FD5TLgc93y54DLZ3eoqsNVdUe3/L/AvcxT\nhVDS5Bs2VE7t6vcA/Bg4daHOSX4XOBv4Vl/z+5PsT3L9XIdPfWMteypNgEVDJcmeJAfneDyvkHr1\nan3MW+8jyW8BXwI+UFVPdc3XAa8AzgIOAx+fb7xlT6XJsOgNhVX15vm2JXkkyelVdTjJ6cCj8/Rb\nSy9Q/qWqbuz73Y/09fks8JWlTF7SyjPs4c/NwDu75XcCX57dIUmAfwLuraq/m7Xt9L7VK4CDQ85H\n0pgNGyp/A7wlyfeAN3frJDkjycyZnD8Argb+aI5Txx9LciDJfuBC4INDzkfSmA31fSpV9QTwpjna\nH6ZXO5mq+gaQecZfPczfl7TyeEWtpKYMFUlNGSqSmjJUJDVlqEhqylCR1JShIqkpQ0VSU4aKpKYM\nFUlNGSqSmjJUJDVlqEhqylCR1NRQX32g48dNd06zffchHj5ylDM2rGPbxZu5/Gy/v1wvZKhoUTfd\nOc21Nx7g6NPPAjB95CjX3ngAwGDRC3j4o0Vt333ouUCZcfTpZ9m++9CYZqSVzFDRoh4+cnRJ7Tq+\nLXvZ067fA9130d6VZO9Sx2u8ztiwbkntOr6NouzpjAu70qZTxzheY7Lt4s2sW7vmeW3r1q5h28Wb\nxzQjrWTLXvZ0mcdrBC4/eyMfffur2bhhHQE2bljHR9/+aj+k1ZzSKyx4jIOTI1W1oVsO8NOZ9Vn9\nfgg8CTwL/ENV7VjK+G77VmArwKZNm1734IMPHvO8JS0syb5ZRxUDW/SUcpI9wGlzbPpQ/0pVVZL5\nEur8qppOcgpwa5L7quq2JYynC6IdAFNTU8eehJKW1UjKnlbVdPfz0SQ7gXOB24CBxk8aLxTT8WwU\nZU/XJzlxZhm4iF+XN110/KSZuVBs+shRil9fKHbTndPjnpo0EqMoe3oq8I0kdwP/A3y1qr620PhJ\n5oViOt6NouzpD4DXLmX8JPNCMR3vvKK2MS8U0/HOUGnMC8V0vPMu5cZmzvJ49kfHK0NlGVx+9kZD\nRMctD38kNWWoSGrKUJHUlKEiqSlDRVJThoqkpgwVSU0ZKpKaMlQkNWWoSGrKUJHUlKEiqSlDRVJT\nhoqkppa97GmSzV2505nHU0k+0G37SJLpvm2XDjMfSeO37GVPq+pQV+70LOB1wC+AnX1dPjGzvap2\nzR4vabKMuuzpm4DvV5XlBaVVathQObWqDnfLP6ZXjmMhVwI3zGp7f5L9Sa6f6/BJ0mRZNFSS7Ely\ncI7HZf39qleUed5ypElOAP4E+Le+5uuAVwBnAYeBjy8wfmuSvUn2PvbYY4tNW9KYjKTsaecS4I6q\neqTvdz+3nOSzwFcWmIe1lKUJsOxlT/tcxaxDny6IZlzBr8uhSppQoyh7OlND+S3AjbPGfyzJgST7\ngQuBDw45H0ljtuxlT7v1nwMvmaPf1cP8fUkrj1fUSmrKUJHUlKEiqSlDRVJThoqkpgwVSU0ZKpKa\nMlQkNWWoSGrKUJHUlKEiqSlDRVJThoqkpgwVSU0ZKpKaMlQkNWWoSGrKUJHUlKEiqalhayn/WZJ7\nkvxfkqkF+r01yaEk9ye5pq990VrMkibLsO9UDgJvB26br0OSNcCn6NX92QJclWRLt3nRWsySJstQ\noVJV91bVoUW6nQvcX1U/qKpfAV+kV4MZll6LWdIKN1SJjgFtBH7Ut/4Q8PpueeBazEm2Alu71V8m\nWY2Fx04CHh/3JJbJat231bpfm4914KKhkmQPcNocmz5UVQtVJFySqqok85Yz7S97mmRvVc37Gc6k\nWq37Bat331bzfh3r2KFqKQ9oGnhZ3/pLuzaApdRiljQBRnFK+dvAmUlenuQE4Ep6NZhhabWYJU2A\nYU8pX5HkIeD3ga8m2d21P1dLuaqeAd4H7AbuBf61qu7pfsWctZgHsGOYea9gq3W/YPXum/s1S6rm\n/RhDkpbMK2olNWWoSGpqIkJl2NsBVqpBb1NI8kCSA0nuGuZU33Jb7PlPzye77fuTnDOOeR6LAfbt\ngiRPdq/RXUk+PI55LlWS65M8Ot91X8f0mlXVin8Av0fvYpz/BKbm6bMG+D7wCuAE4G5gy7jnvsh+\nfQy4plu+Bvjbefo9AJw07vkusi+LPv/ApcAtQIDzgG+Ne94N9+0C4Cvjnusx7NsfAucAB+fZvuTX\nbCLeqdTwtwOsVKvpNoVBnv/LgM9Xz+3Ahu76pJVuEv/bGkhV3Qb8ZIEuS37NJiJUBjTX7QAbxzSX\nQQ16m0IBe5Ls625XWIkGef4n8TWCwef9hu4Q4ZYkrxrN1Jbdkl+zUdz7M5BR3Q4wagvtV/9K1YK3\nKZxfVdNJTgFuTXJf938YrRx3AJuq6mdJLgVuAs4c85zGYsWESi3v7QBjs9B+JRnoNoWqmu5+Pppk\nJ7234ystVAZ5/lfkazSAReddVU/1Le9K8ukkJ1XVpN9suOTXbDUd/ix0O8BKtehtCknWJzlxZhm4\niN732Kw0gzz/NwPv6M4onAc82Xf4t5Itum9JTkuSbvlcev+2nhj5TNtb+ms27k+fB/yE+gp6x3K/\nBB4BdnftZwC7Zn1S/V16n9R/aNzzHmC/XkLvy6m+B+wBXjx7v+idcbi7e9yzkvdrrucfeA/wnm45\n9L6w6/vAAeY5k7cSHwPs2/u61+du4HbgDeOe84D7dQNwGHi6+zf27mFfMy/Tl9TUajr8kbQCGCqS\nmjJUJDVlqEhqylCR1JShIqkpQ0VSU/8PDA5IY+jUjAoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113a6e5c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_train = np.random.uniform(-1, 1, 3)\n",
    "y_train = 0.5 * x_train - 0.3 + np.random.normal(scale=0.1, size=x_train.shape)\n",
    "X_train = PolynomialFeatures(degree=1).fit_transform(x_train[:, None])\n",
    "\n",
    "plt.scatter(x_train, y_train)\n",
    "plt.xlim(-1, 1)\n",
    "plt.ylim(-1, 1)\n",
    "plt.gca().set_aspect(\"equal\", adjustable=\"box\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class AnalyticalBayesianRegressor(object):\n",
    "    \n",
    "    def __init__(self, alpha=1.0, beta=1.0):\n",
    "        self.alpha = alpha\n",
    "        self.beta = beta\n",
    "        \n",
    "    def fit(self, X, y):\n",
    "        self.w_cov = np.linalg.inv(self.alpha * np.eye(X.shape[1]) + self.beta * X.T @ X)\n",
    "        self.w_mean = self.beta * self.w_cov @ X.T @ y\n",
    "        \n",
    "    def predict(self, X):\n",
    "        y_mean = X @ self.w_mean\n",
    "        y_var = 1 / self.beta + np.sum(X @ self.w_cov * X, axis=-1)\n",
    "        return y_mean, y_var"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "analytical_model = AnalyticalBayesianRegressor(alpha=1., beta=100.)\n",
    "analytical_model.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class BayesianRegressor(bn.Network):\n",
    "    \n",
    "    def __init__(self, w=np.zeros(2)):\n",
    "        super().__init__(w=w)\n",
    "\n",
    "    def __call__(self, x, y=None):\n",
    "        self.pw = bn.random.MultivariateGaussian(np.zeros(2), np.eye(2), data=self.w)\n",
    "        self.py = bn.random.MultivariateGaussian((x * self.w).sum(axis=-1), 0.01 * np.eye(len(x)), data=y)\n",
    "        if y is None:\n",
    "            return self.py.mu.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = BayesianRegressor()\n",
    "optimizer = bn.optimizer.Adam(model, 0.1)\n",
    "optimizer.set_decay(0.9, 100)\n",
    "for _ in range(1000):\n",
    "    model.clear()\n",
    "    model(X_train, y_train)\n",
    "    log_posterior = model.log_pdf()\n",
    "    log_posterior.backward()\n",
    "    optimizer.update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sample = bn.sampler.metropolis(\n",
    "    BayesianRegressor(model.w.value),\n",
    "    (X_train, y_train),\n",
    "    1000,\n",
    "    downsample=10,\n",
    "    w=bn.random.Gaussian(np.zeros(2), 0.1)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASMAAAEKCAYAAABZgzPTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8lFe9+PHPd7ZMMtlXAgRCIKwt0DZQWu1Ki7Qqba9V\n22vrcm1763bv9ef1d/vzXper9/e76nW3Wq21akWrol3QtpZCF6iUfStLgEAICUv2TJbJZLbz+2Mm\naUgTMgHCPJN836/XvDLzPM+ZOU8mfHnOOc85XzHGoJRSiWZLdAWUUgo0GCmlLEKDkVLKEjQYKaUs\nQYORUsoSNBgppSwhocFIRB4XkQYR2TvEfhGRH4hIlYjsEZHL++1bLiIHY/seuni1VkqNhkRfGf0S\nWH6W/bcA5bHHA8AjACJiB34U2z8XuFtE5o5qTZVSoyqhwcgYsx5oOcshtwFPmKhNQLaIFAOLgSpj\nzFFjTAD4XexYpVSSciS6AsOYBNT2e10X2zbY9isHewMReYDoVRUej+eK2bNnj05NlVJs3769yRhT\ncC5lrR6Mzpsx5lHgUYCKigqzbdu2BNdIqbFLRGrOtazVg9EJoKTf68mxbc4htiulklSiO7CHsxr4\ncGxUbQngNcacArYC5SIyTURcwF2xY5VSSSqhV0Yi8iRwPZAvInXAl4le9WCM+QnwPHArUAX4gI/F\n9oVE5NPAi4AdeNwYs++in4Dq4/UFWLOrlmULS8hKcyW6OioJJTQYGWPuHma/AT41xL7niQYrZQFr\ndtXy2LpKAN5/9fQE10YlI6v3GakksWRmEXtqmlkysyjRVVFJSoOROi+9zTN/MMyWqkbmT62nJD89\n0dVSScjqHdjK4vqaZ8Zw39LZLFtYMnwhpQahV0bqnHl9AfzBMPdcM4MVi6dpx7U6L3plpEbM6wuw\nauMRVm89xsr1h3G7HMMGot4yXl/gItVSJRsNRmpQvcGjtqmzL4j0BaEt1SNumvU259bsqh32WDU+\naTNNDao3eOypaWZLVWPf9sfWVXLPteV9QSieppnXF8AfCHHPteXap6SGpMFIDao3aCyZWcT8qfVn\nBJGR3ti4ZlctKzdUcd/S2dqvpIakwUgNKivN1XfzYv+h+nO5obE3kOlVkTob7TNS52QkHdK9gU2v\nitTZaDBScRkYfLRDWl1oGozUGYa64ukffIbrkNZhfHUutM9InWGoCa/9+316O6QXzxh8QT+dNKvO\nhQYjdYZ4OpuXLSxh+9EmtlQ1snrrMe69bibw1jy13smy2mGtRkKDkTpD/1E0eCvAtHX18MdN1fgD\nIe69fhbzJmezs7qJfbUteH0BstJcekWkzosGI3VWvQHmsmn5APiDEVZtPML1l0zi0CkvW6oaWbOr\nlvdfPV2H8NV50WCkzqr/zY8bdtdSt7uaVRurqCwv4MFPL2fmxGz8gVDf1ZFeEalzpcFInVVWmoul\nl07ikW+uZs/ja2k51QrAemDLw89R/p5F7C2bgtvl0ECkzkui18BeDnyf6DrWjxljvj5g/+eBD8Ve\nOoA5QIExpkVEjgEdQBgIGWMqLlrFx7CBa1l3tnXxmeu+TMObNeTPnszXfvqPTCqfQHenn6e+9xyv\nPPk6kyrqWfxPS1m18Yiuga3OWcKCUb8U1TcTTcK4VURWG2P29x5jjPkf4H9ix78X+Kwxpn8G2huM\nMU0Xsdpj3uot1azcUIU/EOI9c4p4aPl/0XSgDrn9Hdzw4M0sWTav79hP/PRBgiWFrP/6n3j4879m\nT3kpoJ3X6twk8qbHkaaovht48qLUbIw7602JIgAEuvx89rovU1N5kqyPL0cuK0fEdsZ7fG3VNl5P\nzaDs1gp2rXyVd6XbWDKz6G1LjygVj0Q200aSojoNWA58ut9mA6wVkTDw01jmWBWHsw3Br1hUSopd\n2PHfqzh55DR8eBneidGbG91O21trXgdCvHk82n+U8t4lpGw7TOXKV3hjcTk/f/ngGUuP6JWSikey\ndGC/F/jbgCbaO40xJ0SkEHhJRCqNMesHFhSRB4AHAKZMmXJxamtxgw3B9+8rCryym+0v7OTj3/4I\nTbNL2XK4gcXlhaxYPO2Mof5lCyaz93gLNqeDwKI5VD+3ian+bu5bOnvQpUeUOptENtOGSl09mLsY\n0EQzxpyI/WwAniba7HsbY8yjxpgKY0xFQcHg0xfGm4Gz6L2+AN96dhePratk5W838uv//ANLP3QN\nH/yXd1OUlUq9txunwEPfX8PxndXMSBF2HDpNW1cPJ1t9zCjOYvKN83Glutj81Gbef/V0SvLTdaa+\nGpFEXhn1pagmGoTuAv5+4EEikgVcB9zTb5sHsBljOmLPlwFfvSi1HoPW7KplS1UjFaW57PjmKnKL\nc/jMwx9HRLh5wWR2bjzEHx9+kXBPkIOxMja7jVBuKpdNyyPFYedEZ4CSxTNZ+/u/cc/X7yE3w53Q\nc1LJJ2FXRsaYENE+oBeBA8AfjDH7RORBEXmw36F3AGuMMV39thUBr4vIbmAL8Jwx5q8Xq+5jzZKZ\nRVw2LZ/QtkMc31fLZZ+8lZr2AP/xm8088sg6Nr24m4jbRebl03FcPZfbPnY9BcXZbH5+J9t31+J2\n2rhv6WyKFpXT3dLJyj9sBnT2vhqZRKe3fluKamPMTwa8/iXwywHbjgILRrl6Y87Ae4h6vbrvJDsq\nT2F+sZaMWZNYG7Sz/5md1L5xENPoZfmtCyldNIOrZk9g06FoP9DNS6bzvz73GyKVtSy4/xrmz5xA\nw7EGAE6+WcOqjUfwB8OsXH8Y0E5sNTxdz2gcGbggWu+Viz8Yhk37odPPx//7Hu6/aQ5pTW2YRi+2\n0iL2OVxsPHSaZ7ccix4L/O1IE5FZJYR6gjy+8m94fQEyinNwZaRSs/uYJnZUI5Yso2nqAhg4itYb\nnO6cPxHbG/sIzy7BV5jDnKwUfvLGYSQ/k9RpRZxq9VHb0sbWmhowhmA4zJF6L+Jxk1acw+HdNdz3\n8Mu094TxTMyj6Wg9Sx4o0MSOakT0ymgcGTiKtmxhCfctnY3ZfICIP8idX3gfC6fk8IUvP4U7w032\n/FK6gwEaeo5T76+hLdBAW7CRH7/6Ckcbm7ALvPvWBXT7ArTVNeN22pi/aDqu1nY+t2KBBiI1IhqM\nxrGsNBe3V0xh7S9eZtEtC7nr7xbxle+8SKfXR2BaMS2Bbhp6aglFAiybO4/v3X0HUzOnEgiF2ddw\nmO6Qn+ermnG4HEhbF/5ghHB2OoFOP7Ye7bRWI6PBaJzb+MxWWuu9rPjkcp567SCn9h5HCrOZWJpN\nU08dAKVZ07j/+kWs39dAOJTCtMwyHOKgPdREIGIIZ6QRaulg8YwCrrm6HID6Y41n+1il3kaD0Tj3\n9A+fZ8K0QsqvmcvWV6NzlG1TC9lXX40xhoKUyQSDDv7n2V1ke1LITHXiD4LHkY0/7OPK8hyKS3Kh\nJ8iHrp7O3mYfAHXHdf6yGhkNRuPYwW1H2Pe3g9zy4Lv4+u+3ULnzGDMWTMGWEaA75CPbVQhO4XTJ\nUTbkvcFPOl9gf/5+stLt5KXmALCh6gj1IQPAI09t56+Ho0Fo4/ZjiTotlaQ0GI1ja375Ci63k+65\npex4/RAA73v/FZzoOE2qPZ20LBd1Uw7RldJBRnsuWb5s/O4uDqYfoScEbpuHjp4OFswqBmD/kQZm\nTssDYFJ6SsLOSyUnDUbjVLPXx4u/2cDC5Zfx5sk2IqdbyJ6Sz7df3gwGPGnZ1EysAmBybTmFDSUs\n4RJymyfQmdlGR2YrKfZUukM97DwVnb9sAiGcqdFpIAeqtc9IjYwGo3Fq5a9fp8fro23aRPbuOAbh\nCG25aZzwNpHtzqF5wkkMcFX3AtJD6QCcbOlkccpMim05dE1oIM2ZCkBQgojdhscGE/M9AESCob7P\n0mkhKh4ajMYhry/Aodf2IQ47991/IzmdPtzZHjpTugEhmB2h29NJbkMxDc1dTLr0VcqXvIgto5aq\nU+28I2MmXfjpdgcBmJDjwthtdHb10NQRDTizJuX0fd7qLdU8tq6S1VuqE3G6KknoHdjj0Is7j3Po\n5T1QOoGtB0/TfLqNCVdMpSZUT6ornZaCetw+D5McAcpueBpXahfhkINZRXW0181jsu1OAIKpPgTB\nbo+Qme4Gh41TrZ0ApLr6/WnFVo/s+6nUIDQYjUMF7Z3Q2sni+5chLR0AdOXYoM1gK4oQsYeZ3p3N\nzKueo8eXzv5Xb6O7I5uSSzZTOG0fB/cswCFOAik92MTO0fo28sIRCIY53RRdXMGZ8taf1opFpbid\ndp2jps5Km2ljSLx9M39+4jVw2JizdD5/23CQ4tJ8jnmbyE5NpzPXS1pXOlPKd2AiNg5tvAWfNx8T\ncXDiQAWRiA2fey/OnhQCzh5s2AmZMP5QBGOAUHQibdhm66vLwGkoSg1Gg9EYMnBW/lBadx1l8uXT\nWVCSw/HaZuZfVUrIBEgtdNIjAeZ7IKuojpOVlxP0e/rKhYMp+NrySM+txxlxEnaEQARjDBhDXqab\nD14Rvfo52NSp/URqRLSZNobEk1666UQzdQdP8o/f+jCb3ogO3W9ojA7Dn3K14Ym4mDJ9N9KTR2P1\nXOYVHeWDl71EV7CYb6xbQVdrIQWlleS4Uumwe3HZbYTC0cR3U4syuaYsj98BzvRUiKD9RCpuGozG\nkOHSS3t9AR752SsAlF05k289vgFJT6XW20JWpod6TyNXprqwuds4vPlGHrrxCa4p200wbMdpr2TV\nzivo7s7CZg/THegCl2FCThqd3WHCgXYau3r47u+2AGDLSOWe+aWsWFQ65KJuSvWnzbRxZM2uWtav\n2QNOO7/d30B9XTMpxRmETICMYicAeZ4TBPxpzHG3c03Zbv6w60Y+8uQX6exx8/GrXiPblQuA2xnA\nYGjwdlFWmIlEIjR0BjhypB6AV6pb8QfDrNlV2ze0P1zzUY1vGozGkWULSyjs8uEpLWL3gZOYiMGb\nFt3X6OggLeggPbeOlrpS7lq4lhPefH617VZau7PYeGw+03IOcKIpem+RzR5CjOALBPF2BTARQ2fY\n4OrqBgEy0zhS3xFd8VFEV3xUw0poMBKR5SJyUESqROShQfZfLyJeEdkVe3wp3rLq7dIcQvvR01x7\n4yUsyIvePe3OtiE2G15HB1NT/NjsEaaEhZkFdfxx941EjB2Ao80TyUntJNsdHbo3tghibERMmI7O\n6N3WV80tZoKJkJqXyb03zuZTy+dx39LZrFhUqqNpalgJC0YiYgd+BNwCzAXuFpG5gxy6wRizMPb4\n6gjLjmsDh/oPbz+Kv6sHX3EeDSda8GSn0dzdgT3DAQITM1uJBFNZPnU3Lb4M1h1eBECqy0ZjV/SO\n6mx39L6kEAZb2I7B0ObtiR7ncVNz+BQZxTm4nXYydUhfjUAiO7AXA1WxTB+IyO+A24D9o1x23Ogd\n6vcHw7iddnyv7AXgtc4w4WONhAvSCUb82NOdYAy5eQ3QnMui+Vv59fZ3EcEJGDLdTiQ2KuZ0RgNb\njxHsoehVU7o4iACOtBQyAkE8hRPO+FztuFbxSGQzbRLQv0ezLrZtoKtFZI+IvCAi80ZYFhF5QES2\nici2xsbxNZO8d41rjOGxdZW8/tKbTJhexG1LyiAYJpQd7bQOpQUpdgh2l5+ptujaRG8cu5RwJPq8\nvr0Hu0SbYsYR/ekDbLFgFOiM9iOlelx0NXmZPWfiGZ+rHdcqHlYf2t8BTDHGdIrIrcAzQPlI3sAY\n8yjwKEBFRYW58FW0rt6hfq8vgNvl4JlfPMecxTOYm5PKU4Ary4EJGoLuADOyor+aAlsXgZCDxs6J\nZLgddPijwSfP4wXAZzOkhVz4bGE8ARthIth9IcRh53B1I5FgmI6UlDM+VzuuVTwSeWV0Auj/Vzo5\ntq2PMabdGNMZe/484BSR/HjKqrdkpbl496UTaTjWSNn8Uqqq6rHZBW+4B4/HTcQexhbxEg45Kclo\n5nhbEZ0BIRCb2gEwIaMZXyAFSeuh2+cBG6SEHQiCJ2QwbieV+6JfwaJFZX2fq31GKl6JDEZbgXIR\nmSYiLuAuYHX/A0RkgsQ6K0RkMdH6NsdTVp2p9uBJAEpmT+J4bTOTJuWS5obcvOhiaC6Hj57ODArT\nWzndHl2tsSf01oXktNyTVLcU4/S00OWLrm/U0x7GYXPR3eqjdGo+JSZ6FXXabvULbmVFCfurMcaE\nROTTwItEZxM8bozZJyIPxvb/BLgT+ISIhIBu4C5jotMxByubkBNJEof3HgcguySf6md30u1w0Obz\nUZJXBAZSU7voac8l3+NlR93sM8raJML0vBO8fOQyXHmdeOsngwH8EdLsqZiedkxqCmVdndS6HHSl\npPCzlw4ABkRwO+2sWFSqV0jqrBL6X1is6fX8gG0/6ff8YeDheMuqoW3cdASAr6w5QOtpL5HJORig\nqqUNssHl6qHTn0qKI4A/dOb61TPy6khz9XCkKxfyOmkKpeLqcRMOhynNzOA0jdT1RGh87QCmMJu/\n7Hx7h7XbaT/rVBWl9A7scSI7FETSUqId0hGDIz06kuaXEI6QA7uzh3DAhd1mcDndZ5RdMCm6WH9N\nxEUkYuN0yI4nEL1psu5oGwDFk3LoqW2E4jwunZLDnUvKmDMpG4BLp+RqJ7YalgajMa73xsfW021k\nFGVDINqvM7UkC4CII0IqdkQgHHbiD7rAtFOYGQ1IDluIZbO2cKixBGdhPT0dxfQAtk4bgpDaFiC7\nMJNL0x3QE6T0sml88f0V3H/zHK4oywdgwdRcbaKpYWkwGuN6b3ysqW5gSmkBZbnRyWhBiXZOu1Js\nEBs0M8ZGU1cWBZ42fMEQNoHbL3mNyVmNrKq8GldaBw3eCWDA3xLEJW6CLV1MLiticmd0msjlN1zS\n99krFk+LTgdZPO3inrRKShqMxqjapk6++OQW5pXkct/S2ZguP6cDEY7UtQJwuCk6raMnFMJEojcv\nii1CTWsxswtr6OwOkpvWyt9fvoY3jl3CcWeESNjOoU43Ln8q4VCICSkeMIbiqfnse3UvGRNzeeZQ\nU99Njjq0r0ZCg9EY9ehL+9lS1ch3/rybJTOLaG1spylscJgIADMnZfUdGzF2ImE7DmcPm2rmUZDe\nxndv+z5fvPkX2MTw6Ob3kjv5CM2nS2gWP1Nt0Xlq1xYWIAK3L53L9nVvUnTFDO65tlz7h9Q50WA0\nRj1w81xK8jzUNnfx0xf3Eu4O4EhLIRiItskWlOUjgE2EVJedSDAVp9vHGzWXsnrfOynJbiTf4+Wn\nb9xOMLsNhytATdMEABpr/EzOzaK2spH84hz+/MJuQv4gR1LcYIxeCalzonenjVEl+el8+6NXs2ZX\nLQuKM9kM3Ly4jA0dYbxAKGIwgITtdEuQ1MgEItmt1Es6j2x8H49sfB8ANkeAedc/g689h2NicPs8\n+H1+SgqLqTxQiZQUsm/jweiHTikCkTNWdgR0lUcVF70yGmP6LxvS22fjjESbZiXF2ZTmRRfYP3gy\nOiRvD9sJSpDjJ1Jxelrxh/x4+qUZKpm3GVdaBwcPzSfkDJHdkQGAq10wBiQ3Azl2Gk9RNve+ZwEr\nFpWekRgg3iQBSmkwGmMG+8f/fGzh/f2n2/sWyD/R3A2APegk6AzSemoqYoswae42Up3RDu3s4mMU\nlB7i9OH51NgjuIIuOuv9uGwpVG49DilOFs4ppnl3NTfcvoh7r59FVpqrb7WAJTOL8AfD3HPNDO1H\nUsPSYDTG9AaC/v/4uzqji59FHA7mlEQ7nzu6o/cb5drSQAytXZmcrppHUdl+Mmb9ldLL1jO94mW6\n2vKoOl6O3+3D05xDKBIg1XigtQtbYTZ5Ta34u3q45Kb5b7si23SonpXrD+N2ObSJpoalfUZjzGAZ\nQlyx/3JsDhvHYhlf7bGvPieUTjXQld7Oif2LcKb4ySs5QiRsp7FmNicOzaexqA5H0AlNBkFIawkT\nAfLKJtC55QDpOR5e7oqwbV0lEA2Ia3bVsmRmUd9rpYajwWgcuGZOMc8BLoeNDdVNAJRmZ9DcCKca\nu0lJT8OX0UZeayHVO66jbu9VOO3g606hobCWQIqfwtoptIeaSHdkEznVhmSm0dodYPOftzF9+RVs\nq2nlsml5+AMhVm+pZuWGaNNQ56OpeGkzbQzr7cz2x5YCmTsph/LJ0WZaU6MPQQhG/GR58/CndNPh\n8QJCMJCCrzuFjoxW2rObudJdTqS9J7p2kS8FfD1kTyuCQ3UQDNNWNok7l0TXMFq5oUqzgahzosFo\nDFu99RiPratk1dZjAOyvbuR9180CIOgP4bS5CRo/Ge25pATcnJ54DH9aB2FbiMaCOuon1OD2eTi9\ny0ZnsINCTwHZbT5cbic3LZ1HxrFTkJ5KS3Ymx5s62FndzOIZBZoNRJ0TbaaNZSZ6RTRtcg77gMkZ\nKaw5cBqACW4HeDI52VFPxISZ2TWdSlNF3eQjsbJQ6Csk7WQhzYGT2LCTHfbQdrwWW0kBT68/hOyp\nZtZ7KlhwdRkgzCzOYsXiaRqE1DnRK6MxrHei6gdunAPAS1ur2XWqHbHbaGrwEg5GZ+Ybm4/7F13G\n4vaFFJ+YRk5zIfObL+FK5tLZ00Iw4ifbVYj34CmwCXnlEyiqOUkkGOL+/72CbE8Kf9x0lEOnoutk\nD0yRpFQ89MpoDOsdWYtEIohNOHWyDdtcgbQUAm1dpJYVk5GSRpOvmV+8coCO7jAesvB0ZeED6rOa\n6Qy1MSmjANNmJ9zoxVZSQGtPBLNmB0zO51DYxrKFJeypaWZLVWPf/U2PxUbWtANbxUuD0Thgs9nI\nKsgix+MgbVI2tfU5tB46RTAYZlZBCdvqDnK0tYa5hdPJz0glPzOFNfsOsOdUE+lODyaYReRINTgd\n2EoKoPI4tHRw3T+/t2+ax7/etvCMKSCgQ/pqZBIajERkOfB9outYP2aM+fqA/R8C/o1o9vYO4BPG\nmN2xfcdi28JAyBhTcRGrnlS8vgD2rDQ8gSDpqU7anS4whvI0O3det4DGv/io8day8/SbuBs9RAgS\nCAdIs2eS6SiABi+m3cfkxTPImZrL/l/+lfTiHP7PF1Zgt0fv1h54f5NeEamRSlgw6pei+maiSRi3\nishqY0z/rLDVwHXGmFYRuYVo/rMr++2/wRjTdNEqnaRWb6mm2eGkY18toapGJpUW0FR9ikmRME+8\neohwMJVcVzFud5DOHh/+oJ08VzGpjgyMr4dQ1UkkM43TKW5mt3Wwv7aRD377I3T2hFm9JTqUrwvu\nq/Nl6fTWxpiN/Y7fRDQ/mopT7+x5b3cQ8rIIVNYioTCn2v3MnDuZV187wOSl8wFIc2Rw56Iyli8s\n4Uu/28rJVh8OE6F7fw12hx2ZHW1y7f/Nq2QXZHLjvdfxrWd3saUqmqVXF9xX5ysZ0lv3+jjwQr/X\nBlgrIttF5IGhCo3n9Na9k2a3HG5AJuZCJAKnW8j1uPDmZiJ2OzVvHMLEbgHoCYT4zz9s42Srjwyn\nDf+bNdDdw3V3LOLyORPh8Anqdx7l7i/8HS/uPcmWqkaKslK5c0kZS2YW6QiaOi9JMbQvIjcQDUb/\n1m/zO40xC4FbgE+JyLWDlTXGPGqMqTDGVBQUFFyE2lrHsoUlLJ5RQL23GybFzr22kZ5QhPruEEyb\ngGn3kXK8gXcvnMyxxk5qm7twBwP4tx3GeDu5/Kb5TJ0xgQdumEnW67spnjGB935iWd/s/3pvN9ke\nF5sO1etSIeq8JLKZFleKahGZDzwG3GKMae7dbow5EfvZICJPE232rR/VGieZ3lGu1VuPYSIR1jyz\nHmdLK5fMnsCa3XWYgizc3X66ahp47mfryJiYS6ixnY4OHzgd3P4PN5BTnMPK9YdZ+8PnaKtp5PPP\nfQGny8mKRaXRDzFGR9DUBWH19NZTgKeAe40xh/pt94hIRu9zYBmw96LV3OJ6bzqsbepkza5aViwq\n5cM3zObaOxZzeudR3CZCWVE6GalOPvngUoqWzCKUmkLL0XqcdhszK6bzwU8uY8kV03h2azWmycvJ\n1ZsoWzqfWdfNY9XG6F3a9143s28NI118X50vq6e3/hKQB/xYos2C3iH8IuDp2DYH8FtjzF8TcBqW\n1NtX1HsjIkSH2q/7wNWs+vafefbnLyOLoimsNx9uoMXl4rJbL6fB66Pe6+cd15Zz73Uzue/Hr9Le\n7sf21AZIcfGZ7/1D33v3vqdSF4rV01vfB9w3SLmjwIJRr2AS8voC+AMh7rm2nOvnTWT+1Pq+ptPM\niunkzZxI88Z9XHr7EhaW5XPF9EIAsj0pvHm8BYB9tS14fQHuv2kOX/377xE42YTcvZQDXn/fe2lz\nTF1oSdGBreK3ZlctKzdU4XbaKclPP6PpJCLc/193Q0sH82pPsmLxNJ58/TBbqhrZfSzaHedJcbCz\nupn/fmoHf/zG0wR2VLHwo0u5/5M3991trc0xNRp0OsgYM9yVy9I7l7Dpg1ez6v/9iV3OVA6GhUun\n5NITDFPv7WZibhqH6trY/sPnYGslM5Yt5N9/8DGyPSkX8zTUOKRXRmNM74L4a3bVDnnPz2cevo/U\ngiwO/MevMHurMSbSN+O+xO/Hs+oV2FrJvA9ew8PPPaSBSF0UemU0Bg3XyZyZl8E3XvoSX7zjmzSv\neo3aXYeJYCMtFGRtdT0pHjfT//EWPve1D/TNPRuof240bbKpC0GD0RgUTydz+ZxJ/Gbvd3j2R3/l\n5d//DXuLj5z8dK755LvYk53NjpPtbDpUT0l++qDldVRNXWjSOxVgPKioqDDbtm1LdDUsbdXGIzy2\nrpKSPA9f/kDFGcFIM8Wq4YjI9nNdQUP7jJLUua6mOFi52qZOvvjkFmqbOvumkNQ2d7HpUP0ZZfsn\niNRRNXWhaTMtSZ1rM2mwco++tD92c+R+vnb34kEXSoP4mn9KnSsNRklqsMAQT6fyYOUeuHkusD/2\nc/BEkGfbrtSFoM20JDVYM6l/M2owA4NVb5MtM83F1+5ePGRntVIXg14ZjSHDNaMGNtH6v+69N0k7\npFWiaDAaQ/o3o2qbOnn0pWjTq/eKZ9nCEvyBEP5gGK8vcEbw0qF6lWgajMaogZ3SEA1WbpeDx9ZV\n9i0T2xuslI3dAAALeUlEQVR4tHNaJZoGozFqYKd0r6GCjnZOq0TTmx6VTu1QF4ze9KjOMNIbIocb\nhVPqYtBm2hg00s5o7S9SVnDOwUhEvmSM+eqFrIy6MEYaXLS/SFnB+TTT3rYc7EiJyHIROSgiVSLy\n0CD7RUR+ENu/R0Quj7fseKbzxlQyOuuVkYi0D7ULSD2fD44zvfUtQHnscSXwCHBlnGWVUklkuCuj\nNqDcGJM54JEBnDrPz+5Lb22MCQC96a37uw14wkRtArJFpDjOskqpJDJcMHoCmDrEvt+e52fHk956\nqGPiTo09ntNbK5VMzhqMjDH/YYzZIiIrReR+EZndb9+/na2sVYzn9NZKJZN4O7B/DhQDPxSRoyLy\nJxH55/P87HjSWw91TFypsZVSySOuYGSMeQX4v8AXgZ8BFcAnzvOzh01vHXv94dio2hLAa4w5FWdZ\npVQSies+IxFZB3iAN4ANwCJjTMP5fHCc6a2fB24FqgAf8LGzlT2f+iilEivemx73AFcAlwBeoE1E\n3jDGdJ/Ph8eR3toAn4q3rFIqecUVjIwxnwUQkQzgo8AvgAmAZvdTSl0QcfUZicinReT3wE6i9/M8\nTvSGRJUEzjWTiFIXU7zNNDfwHWC7MSY0ivVRo0BXcVTJIN5m2rdGuyJq9OisfJUMdAmRcUBn5atk\noIurKaUsQYORUsoSNBgppSxBg5FSyhI0GCmlLEGDkVLKEjQYKaUsQYORUsoSNBgppSxBg5FSyhI0\nGCmlLEGDkVLKEjQYKaUsISHBSERyReQlETkc+5kzyDElIvKKiOwXkX39s5GIyFdE5ISI7Io9br24\nZ6CUutASdWX0ELDOGFMOrIu9HigEfM4YMxdYAnxKROb22/9dY8zC2EPXwlYqySUqGN0G/Cr2/FfA\n7QMPMMacMsbsiD3vAA4wRNZYpVTyS1QwKorlPwM4DRSd7WARKQUuAzb32/wZEdkjIo8P1szrV1bT\nWyuVBEYtGInIWhHZO8jjtv7HxdIRmbO8TzrwJ+BfjDHtsc2PAGXAQuAU8O2hymt6a6WSw6gtO2uM\nuWmofSJSLyLFxphTIlIMDJoQUkScRAPRb4wxT/V77/p+x/wM+MuFq7lSKhES1UxbDXwk9vwjwLMD\nDxARAX4OHDDGfGfAvuJ+L+8A9o5SPZVSF0migtHXgZtF5DBwU+w1IjJRRHpHxt4B3AvcOMgQ/jdF\n5E0R2QPcAHz2ItdfKXWBJSQ7iDGmGVg6yPaTwK2x568DMkT5e0e1gkqpi07vwFZKWYIGI6WUJWgw\nUkpZggYjpZQlaDBSSlmCBiOllCVoMFJKWYIGI6WUJWgwUkpZggYjpZQlaDBSSlmCBiOllCVoMFJK\nWYIGI6WUJWgwUkpZggYjpZQlaDBSSlmCBiOllCVYNr117LhjsbWud4nItpGWV0olDyunt+51QyyF\ndcU5lldKJQHLprce5fJKKYuxenprA6wVke0i8sA5lNf01koliVFLVSQia4EJg+z69/4vjDFGRIZK\nb/1OY8wJESkEXhKRSmPM+hGUxxjzKPAoQEVFxZDHKaUSy9LprY0xJ2I/G0TkaWAxsB6Iq7xSKnlY\nOb21R0Qyep8Dy3grjfWw5ZVSycXK6a2LgNdFZDewBXjOGPPXs5VXSiUvK6e3PgosGEl5pVTy0juw\nlVKWoMFIKWUJGoyUUpagwUgpZQkajJRSlqDBSCllCRqMlFKWoMFIKWUJGoyUUpagwUgpZQkajJRS\nlqDBSCllCRqMlFKWoMFIKWUJGoyUUpagwUgpZQkajJRSlqDBSCllCZZNby0is2JprXsf7SLyL7F9\nXxGRE/323Xrxz0IpdSFZNr21MeZgLK31QuAKwAc83e+Q7/buN8Y8P7C8Uiq5JEt666XAEWNMzajW\nSimVMFZPb93rLuDJAds+IyJ7ROTxwZp5SqnkMmrBSETWisjeQR639T/OGGOAIdNOi4gLWAGs6rf5\nEaAMWAicAr59lvIPiMg2EdnW2Nh4PqeklBpFlk5vHXMLsMMYU9/vvfuei8jPgL+cpR6PAo8CVFRU\nDBn0lFKJZdn01v3czYAmWiyA9bqDt9JeK6WSlJXTWyMiHuBm4KkB5b8pIm+KyB7gBuCzF6faSqnR\nYtn01rHXXUDeIMfdO6oVVEpddHoHtlLKEjQYKaUsQYORUsoSNBgppSxBg5FSyhI0GCmlLEGDkVLK\nEjQYKaUsQYORUsoSNBgppSxBg5FSyhI0GCmlLEGDkVLKEjQYKaUsQYORUsoSNBgppSxBg5FSyhI0\nGCmlLEGDkVLKEhISjETk/SKyT0QiIlJxluOWi8hBEakSkYf6bc8VkZdE5HDspyZxVCrJJerKaC/w\nd8D6oQ4QETvwI6J50+YCd4vI3Njuh4B1xphyYF3stVIqiSUkGBljDhhjDg5z2GKgyhhz1BgTAH4H\n9GajvQ34Vez5r4DbR6emSqmLJSGpiuI0Cajt97oOuDL2vMgYcyr2/DRQNNSbiMgDwAOxlz0iMhYT\nPuYDTYmuxCgZq+c2Vs9r1rkWHLVgJCJrgQmD7Pp3Y8zZMsiOiDHGiMiQaav7p7cWkW3GmCH7qJLV\nWD0vGLvnNpbP61zLjlowMsbcdJ5vcQIo6fd6cmwbQL2IFBtjTsVSXTec52cppRLMykP7W4FyEZkm\nIi7gLmB1bN9q4COx5x8BLtiVllIqMRI1tH+HiNQBVwHPiciLse0TReR5AGNMCPg08CJwAPiDMWZf\n7C2+DtwsIoeBm2Kv4/HoBTwNKxmr5wVj99z0vAYQY4bsblFKqYvGys00pdQ4osFIKWUJYzoYne+0\nE6uKdzqMiBwTkTdFZNf5DLmOtuF+/xL1g9j+PSJyeSLqeS7iOLfrRcQb+452iciXElHPkRKRx0Wk\nYaj79s7pOzPGjNkHMIfoTVivAhVDHGMHjgBlgAvYDcxNdN2HOa9vAg/Fnj8EfGOI444B+Ymu7zDn\nMuzvH7gVeAEQYAmwOdH1voDndj3wl0TX9RzO7VrgcmDvEPtH/J2N6Ssjc/7TTqxqLE2Hief3fxvw\nhInaBGTH7i+zumT824qLMWY90HKWQ0b8nY3pYBSnwaadTEpQXeIV73QYA6wVke2xaTFWFM/vPxm/\nI4i/3lfHmjIviMi8i1O1UTfi78zKc9PicrGmnVxsZzuv/i+MOet0mHcaY06ISCHwkohUxv5HU9ax\nA5hijOkUkVuBZ4DyBNcpIZI+GJnRnXaSMGc7LxGJazqMMeZE7GeDiDxNtNlgtWAUz+/fkt9RHIat\ntzGmvd/z50XkxyKSb4xJ9km0I/7OtJl29mknVjXsdBgR8YhIRu9zYBnRdaSsJp7f/2rgw7ERmiWA\nt18z1cqGPTcRmSAiEnu+mOi/yeaLXtMLb+TfWaJ75Ue5x/8Oom3VHqAeeDG2fSLw/ICe/0NERz7+\nPdH1juO88oguKncYWAvkDjwvoiM4u2OPfVY+r8F+/8CDwIOx50J0ob0jwJsMMTJqxUcc5/bp2Pez\nG9gEXJ3oOsd5Xk8Cp4Bg7N/Yx8/3O9PpIEopS9BmmlLKEjQYKaUsQYORUsoSNBgppSxBg5FSyhI0\nGCnLE5ErYqsPVMVmgkui66QuPA1GKhk8AtxPdJpEObA8sdVRo0GDkUo4Efm8iPxT7Pl3ReTl2PMb\nRWQdkGmM2WSiN8U9QXKvUqCGoMFIWcEG4JrY8wogXUScsW1riN7h2ytZZuyrEdJgpKxgO3CFiGQS\nnbrzBtGgdE3suRoHkn7Wvkp+xpigiFQDHwU2AnuAG4AZROffTe53eLLM2FcjpFdGyio2AP9KdImT\nDUQnXe400Zne7SKyJDaK9mE0aeeYpMFIWcUGoBh4wxhTD/hj2wA+CTwGVBGdBf5CQmqoRpXO2ldK\nWYJeGSmlLEGDkVLKEjQYKaUsQYORUsoSNBgppSxBg5FSyhI0GCmlLOH/AxN/MpooF1ubAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113b48160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "w0, w1 = np.meshgrid(np.linspace(-1, 1, 100), np.linspace(-1, 1, 100))\n",
    "w_grid = np.array([w0, w1]).transpose(1, 2, 0)\n",
    "plt.contour(w0, w1, multivariate_normal.pdf(w_grid, mean=analytical_model.w_mean, cov=analytical_model.w_cov))\n",
    "w_sample = np.asarray(sample[\"w\"])\n",
    "plt.scatter(w_sample[:, 0], w_sample[:, 1], c=\"steelblue\", s=1)\n",
    "plt.xlabel(\"w0\")\n",
    "plt.ylabel(\"w1\")\n",
    "plt.xlim(-1, 1)\n",
    "plt.ylim(-1, 1)\n",
    "plt.gca().set_aspect(\"equal\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sample = bn.sampler.hmc(\n",
    "    BayesianRegressor(model.w.value),\n",
    "    (X_train, y_train),\n",
    "    sample_size=1000,\n",
    "    step_size=1e-2,\n",
    "    n_step=10\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASMAAAEKCAYAAABZgzPTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8nVWd+PHP92652fcmTZu26RK6QEmhlAItW2kFRllU\nEEZw+YkdBZxRf85L/DnjODrOOPNzGccF7SCC9ieMqEhVlkKltiylCw3d26TN1uz7dnNzt/P747mJ\ntyFp0iW9T5Lv+/XKK899lnvPkwvfnnOec85XjDEopVS8OeJdAKWUAg1GSimb0GCklLIFDUZKKVvQ\nYKSUsgUNRkopW4hrMBKRx0WkSUQOjHBcROS/RKRcRPaJyGUxx24WkaPRY49cuFIrpcZDvGtGTwA3\nn+b4LcCC6M964FEAEXECP4weXwzcKyKLx7WkSqlxFddgZIzZBrSd5pTbgZ8byw4gQ0SmAyuAcmPM\nCWNMAHg6eq5SaoJyxbsAo5gB1MS8PhndN9z+K4d7AxFZj1WrIjk5+fKFCxeOT0mVUuzZs6fFGJN7\nNtfaPRidM2PMBmADwPLly83u3bvjXCKlJi8RqTrba+0ejGqBwpjXM6P73CPsV0pNUPHuwB7NJuAj\n0adqK4FOY0w9sAtYICJFIuIB7omeq5SaoOJaMxKRp4DrgRwROQn8E1atB2PMj4HngVuBcsAHfDx6\nLCQiDwMvAU7gcWPMwQt+A0qp8yauwcgYc+8oxw3w0AjHnscKVkqpScDuzTSl1BShwUgpZQsajJRS\ntqDBSCllCxqM1Ig6fQGeeeM4nb5AvIuipgANRmpEm0treGzLETaX1ox+slLnyO4jsNUF1ukLsLm0\nhnUlhawrsQa5D/xWajxpMFKnGKgNAdx19TzuunpenEukpgoNRuoUWhtS8aLBSJ0iPcmjtSEVF9qB\nrZSyBQ1G6ozEPu7XR//qfNJmmjojsR3cwCmd3UqdCw1G6owM18Gtnd3qfNBmmjojAx3c6Uke0pM8\nrCspZHNpjTbV1DnTmpECTh3smJ7kGfN1A802fzCM1+1kZXEeO441nvH7KKXBSAHvHuwYayBQXbsg\nl9p9VbTWtRHsD7H6A1cONtH8gRCPbTnCvqpWdpY3D/s+Sp2OBiMFwMriPPZVtbKyOO9dx17cU8VP\n//N5Nm7fh6+la3D/o5/7GX+1fi33/9NdhFwuvB4XK4vzWDq78V39SGdb81JTR7zXwL4Z+B7WOtaP\nGWO+OeT43wMfjr50AYuAXGNMm4hUAt1AGAgZY5ZfsIJPIjUtPWx4+RCzclLYWd7M0tmNFOakDB7v\n6ejltUeexLx5lILl8/jo4w8yY0E+fT1+fvuff+TZ7/2RyoM1fOOPXxqsCcVeP2DTrko2bivDHwxz\n/3XFF+z+1MQRtw7ssaSoNsb8X2NMiTGmBPgS8GdjTGwG2huixzUQnaUNLx9iZ3kzxxu7eWDNwlNq\nNO2NHXzuuq9wbFc5Dz66nh/u+FdWvvdy0gpzeacfPv2TT/HZn/wNeza/wxP/8DRwmmVHjDn1t1JD\nxLNmNJiiGkBEBlJUHxrh/HuBpy5Q2aaM9WsXA4dYv3bxKTWaurp2/u6aL9Pd0IG590beSk7lRn+I\n9CTPqf1Ln1jD0V3Hefrff8fS65dwJCGRjdvL8QdC3H/9RYPvd9uKIrwelw4DUCOK56P9kVJXv4uI\nJAE3A7+J2W2AV0RkTzSFtToLhTkpfP3eFacEokgkwj/e/R06alqY89k7kPkz2FvROriu0ZLCLPLS\nE3n9aAM1LT08+L2Pk1c0je9+/kn6gmHrTURO+ZzYIQFKDWeijDN6H/D6kCbaqmjz7RbgIRG5drgL\nRWS9iOwWkd3Nzc0XoqwTUmzz6pff+C3Vbxzhmr+7jcvXXAJAXnoi1S09PLJxB4//6TCNnX0cPtnB\nV57exf/sqGDO7StpPnyS1kPVrJify/VLCoZ9b6VGEs9m2kipq4dzD0OaaMaY2ujvJhF5FqvZt23o\nhcaYDcAGgOXLl2uHxQgGml5Nh2rY9M+/ouimS5n73uXccMlMjjd08vbxFurLGzA+P1n5GSyYkUV3\nKEJdu4+N28q444aluH78ArVbD1B22SKKC+qs/iERMIaN28sBfdyvRhbPYDSYohorCN0D/PXQk0Qk\nHbgOuC9mXzLgMMZ0R7fXAV+7IKWepNaVFBIJhXnx04/iSk+m8oolVL12HK/HhaOxndDOIxAIAdBS\n1USb08GDn1nH7yvaqWv3kZLq5ar3Xs7+bYf568+8l3cqW9hf3Q7AB1fOZcX83GGHDSg1IG7BaKQU\n1SLyqejxH0dPvRPYbIzpjbk8D3hWrH4JF/BLY8yLF670k8PAY32rExs2P/EqJw+fRO65gcsWFbCo\nIJ2aXeXsePEdCubk8rH7rqGqP8zTLx0gq7WTR7+/GS4p4vJLC7n+4hlUX1RIx6930FnZxP7abgCW\nFeXgdTuGHTagVKx4p7d+V4rqmCA08PoJ4Ikh+04Al45z8SatgQGIe060sLeiBTiEv9dP9a9fg1nT\nWPHey/nkTYv553/5HZWHa1ly5Xy+9sj7yExJoNMXIC3Fy8q5OTz0mSfpOVRF9ooifvjiQd7utjqv\n6/ZVQHYOy4qy+dL7lwHokzQ1qonSga3Oo4H+oXl5qayYn8u9qxZQ9/xu6PGTdefV/M26JfzkZ9uo\nPFyLY04exxISeWXfycHr/YEQL+6vo3DVIgiGeW3LASuoZaVCYgKdxxt4YM1CvvT+ywYn1OqTNDUa\nnQ4yBQ1M/bh52SwKc1L4xQv7aHphNywspCMrg/95cR9vbTmAa1o6FOaS6BaON7ay4U/tbNlfQ1un\nC4dY/44lF2TSXdVM4aJCGrr8uAuyoKljcDa/Tv9QY6U1oylo68E6dpY3s/VgHQB92w9AIMR7Pncb\nH7pyDluf3YXxuGD+DCImRGV3Jb/Y+Trfe3k7BxoqafJX4fGEyEtP5N4ProBQmLUz0ynMTiaQlUb1\n4ZO8tLdac66pM6LBaCqKmZrR0uHjDxte5rJ1l/KFT15PX0Ujvq4+nAsLCTqCNPXXEIoEyHBPY05q\nETkJMzEYTnScoKatnT3NPpweF1vfOEZNay+zi6cT6PHT1drNfdcu0H4iNWbaTJuCYqdmfPtfn6Ov\nrYfkVRfT1NTFs8/uRqZlEEn10OyvwiEOLskrpq3bEAqD1+kmzzubhr5KesIttPRkEklNoqq8kbsf\nLqFJAlQBv3phP1deuzDet6omEK0ZTTFDl/LY98xrkJlCb2Euv3z6TYwxOGZPoy3QABhyE2bS7fvL\n9YluIdnjJtWdiS/ko6atA29WMmFfP0dr2vhzdQcAxWkedpY3azNNjZnWjKaY2EmuSz3QW1ZH5gdW\ncdPC6fz7V7Yz5+JZHHL1Egj2kenJA7dQlV9OX2IvYgRPwEtB7VxSPel0BlvoDXVy3dIi/lxWT5YY\nSPICsGpWJtctm6/NNDVmWjOaIgbmhy0pzBocDb35iVdxelx0FM/iuxteJRI2zL5sJt3BVhKdKeRM\nS6JuThl9iT3MjuSR3JOO39tL07QajHGQlpCKyx0gO9cayNjT6SM/Lx2AIyea9XG+OiMajKaIgRrR\nU6+VsbO8mTcO1fHq069TeM0ipuek4D/ZgkzL4LWaKtwuB1cUF3Ek5wiRiOEqXwl3Zq1gXudcslrz\n6UnroDutnYvy8+jo8/HsvmoAdh6qpyE6Gba7y6cTY9UZ0WbaFDHQXBpYFrags4vuth568nKY3twJ\n4QhpC3I51lnN3OzplKdWEPYbCmvm0yXC7+uqAMj059Gb3EVbdj01J6wmWVCCuJwOHKEwOJ0A7K9o\nYdPOilPWNFLqdLRmNEUMjIJOizab9rxUisvjYuYV86k5fBJ3ehK14U4Eocndx4lAEzM7C/EK5C/Z\nyvwrXyI1pw5ByGjPJeQO4k+0Js46JAxOB+FgiARX9D8pwylrGukyImo0WjOaYjaX1vDfrxzG89Tr\nTCsp4mRTN/T4Cc7LozfUSWZyOs15TWSF08l3+Jl5w0t4EnsJh1xk5J+kpXIRx/evACCQ6EMQPB4I\nOZ0QinDZ3BzeAOZPTwNj6PQF3r06pC4jooahwWgK6fQF8AfDFIcDHG3uZObH1iCRENVAX5YDjMGX\nESAkIbJbk5mz6nn6fSkc2no7/p4MZi55i2lFh2msuAhn0E0woR+HOAlHQngTXHi8Lk7UWY/2m30B\nNm4vx+txcdfV84bNRKtULG2mTRGdvgDfeq6UjdvKcByuwuF20T+vgOrDdSRPSyPi6aMoJwtfVicp\nvlTmL9wLxkHbvjsxvjwiYRe1h5cTiTjImnGCxJCXgLsfB06CkRD+YJhuX4DGlh4A1pTMPmWBf50s\nq0ajwWiK2Fxaw87yZgoyk6h56xgpF81g34kW6OunK9NLT38fSfkJBJ1BZjuCpOedpPbwZdS3ODHR\n6SPhYAK+jmxSshqRoIuwKwQiBEJhMIaEBBf0W31Cc2dna/BRZ0SD0RSxrqSQFfNzqa1qpqemhe4Z\n0ygwVge0P806p97TTnLEw8KLDtDXnUFTxWKW5J3g/9z0KJ++2sqF0Ns+jeTMFhxhJ2FnCAECoQhE\nDFmpXvISrKdpyelJ8bhNNYFpMJoi0pM8rF+7mDn+PgBuvXM57u4+nGlJ+OnDneChLtJGoQRwJ3VS\ne+gyvnTjz/nWbd+npKCM25a8xtysWvy9aTicYbyuIIjBYHA6BI9AVyBMY5211GxmdPCjUmOlwWiK\n6PQF+OGLB6jYcxxngpucufmcKG8klJlIyATILEwAYFZWIwF/Eou8Xaye+w6/Kr2RB371FXr6vdxd\n8gpBv1XjSXAHMBgwhkS3i4A/iC9sSAtbqz260pLjdq9qYtJgNEVsLq1hb0Ur1LfhLMjmyef3W31B\n2W4AWpw9JAVdZOTWkey/mP+1YivNvdN4cvetNPWk8UblUpZOLycSsh7AOpwhHAhOF/j6rT4jXE6S\n/H4QeLu595TP13FGajRxDUYicrOIHBWRchF5ZJjj14tIp4iURn++MtZr1anWlRRy78o5OBraCORm\nkBIMAtDtDCIOB31eH3MT+3E4I6R09DMzvZJf7rmeiLH6gE60FpCZ1EN6QjTIOCIQcRAMh8hOtEZi\np6V6kY4e3JmpXHPxqfk4B8YZ6Sx+NZK4jTMSESfwQ2AtVjbZXSKyyRgzNL31dmPMe8/yWoVVK9m0\ns4K2snoi/UEWXFXM0fouTKKHfvrwZnrokx6mpbQR7PeyOn8Pbb5UtpRdMfgezb2ZAGQmWlk/ghgc\nYRfhSIQsbwKNQHfY0H2iCVKT2HHs1EwgOs5IjSaeNaMVQLkx5oQxJgA8Ddx+Aa6dcjaX1rBxezkv\nPrcbgC/+7Xtw9PiIZCYSMWHc6VaixcysJpztWVwx6wjPH76aYOQv/1YZY03tcLqtZla/EZwhq9bU\n1uoHYG5hJs5uHzOLpr0r6Og4IzWaeAajGUBsnf1kdN9QV4vIPhF5QUSWnOG1kzq99Wj9MLHLhuSl\nJ2JONpNSkIXD7SLcHyKUYfUX9bj95IrBneBnhokA8EblJYPv43EJbqc1DEBc1m8f4Axbwaivy/r8\nE+1+wp29BJO9GnTUGbP7dJC3gVnGmB4RuRX4HbDgTN5gMqe3Hmm+18Bqjh29/fx6RwXLirJp7OyD\nhjayLy2irKwBgLTcBDp6u/F7/MzzWgFluttHIOSiuj1/8P2CIUN2cicAPQ5DSshDnzNMUr8DiOAN\nGEJOB55gkEA4woy5mjlWnbl41oxqgdi6/MzovkHGmC5jTE90+3nALSI5Y7l2KlhXUnjKlIsBA0Hq\neKPVvzMvL50PXV4I7T1kzctjx94qEKGxz4cnwUPEGSbV7SMcclOY2kp1Rx7haMc1WBPw81Nb8QUS\niCT00edLxghIv+BAKHC7cCYncPdiKwjdcuPiC/Y3UJNHPGtGu4AFIlKEFUjuAf469gQRyQcajTFG\nRFZgBc9WoGO0a6eCgX6YoWLXLtpxrJF1JYWUv32Cp4DSnjCevdWQ6CFkgkhCNP9Zgp/+nlSm5dYO\n1ooEKxABFGXVUdE2ndSMTtq7sgAwPnA6PNRUt0F6MvvfrgBgZ3s/NVuPctuKIm2uqTGLWzAyxoRE\n5GHgJcAJPG6MOSgin4oe/zHwQeDTIhIC+oB7jDVRathr43IjNhQbpApzUqhp6eFbP90GQP68PPoq\nGgmkJBA2/bgTrGCRmNhLf1cWOcmdvH1yIbNzU6hqtia9OiTCvOxa/nR8Gc7sHtrqZ1hRyh/BbRIw\n/V04kr2U76sGt5NXqjqQGqtZp4urqbGKa59RtOn1/JB9P47Z/gHwg7Feq96d/QNgw8uHaD3ZCkA9\nDiKtPYRnZAD9GI8BAx5PPz3+RBJcAXzBBJwOweUQQhHD/OyTJHn6Od6bBdk9NPV78fR7CYfDLMhM\np4FWJCWRxLYuTGEuaZmJNHX6T1lcTanR6AjsSWagv+hbz5UOPmVbv3YxuSZMcmYy911bjAlHmFdk\nNbUSUxykOBJwuvsJBzw4HYZwxMmJxm5CEauRdumMYwBURTxEIg4awk7cPivQtdVagyA/dNMieisb\nSZ2bT1OnnxXzc7l+SQG/2HqUX/z5mI68VqPSYDTJDMzOj81ZVpiTwtxkF3mzcrh6bjYAHWFrBHZS\nspM0lwcRCIfd+IMekhOsybRpiW5cjhDrLtrJseZCXLmNdLZNI+QQUkMJCEKgrhuSvewrraavu4+m\n1BRWzM/lC7eXsONYIxu3l7NxW5mOvFajsvujfXWG0pM8fOH2ksGm2oCWunYCCQnUNnQBUN/jBxd0\nBQKkeK3mlDEOWn0ZFKRZ/T1dfUE+uPTPzExv5huvfYiE4oOU1RaBgf72EB7xQmcnnlk5HHndGmJQ\nsKyI9WsXk57kYV1JIf6AteaRjrxWo9Ga0SQ0dLRzpy9ATXUrJ/0htr5j1VByMq35ZL5gkI5ua6a9\n22WoaMtnXnYlYMhNbufDl23mzcqLqUswRMJOKv1JePyJ9PX1kyFeMIbkaRmY47VIVir1xsnWA7WD\nfVe3rSji/uuK9amaGpUGoylgc2kN/V0+8gsyCfqtvhvj+EvncjjiIBJ2gtPPwYYSclM6+M5t3+Mf\n1v4MEcOGt95H2vRy2htn0eMOktSTCECRMwkEuhxOqGjAzI8OghfRibHqjGkzbQpYc8kMNvQHWb1s\nNp68NN4Emnus+WQiggGC/V7cXh8vH7iS6WmruGHe2wTCLn7y5h2EMjpweQI0d88EfCT5XEQcHioP\nNuBOS8K0dmGCIVbeeDHFq+dz2xVzBj9bm2dqrDQYTQEJ0flmSWmJ7K9uAyA/M4nulnYk7MThjhDs\nzSIxrZ2+YCKPvvEBHn3jAwA4XAGWXP87fJ2ZHOgP4A0m09beQ4ojA9PVyLwrF1C+9wQAn/n0GqYV\n5gx+rqYkUmdCm2mTUKcvwC/+fIxfbD1Kpy9AX7QWdLzVx57jTSBw7eICwJrsGiBId0cG3pQOcARJ\n9bqYmWWt6Fi45C08Sd0cLVtKyB0ipcNaFsRrzTThr25aQnpTG0l5GSRkp134m1WThgajSWhzaQ0b\nt5WxcXs5m0tr8Pf2A3DZwnwumZUFBhxiffVpjkSC7iCdDXNwOCPMWLyb+dPTSfS4yJheSe6cYzSW\nL6XSGcETTMDf0k+iy0tSSx8kuNl6opmO/ZX4CvO0f0idE22mTULrSgrxB62lYNeVFNJ8zJpDnJGR\njLvemuLx5tEmAIpSMqmjntaeVBrKl5A//yCtJ/04cTG3oIzejmzKqhbQP6OSyyMLOBipJdlk0d9c\nT0ZxAd7aVgiEmLd6sfYPqXOiwWgSSk/ycP91xYOvG8NWn9H2Iw1YyYXA47C++vqqIMyG3pQuag9d\ngTvBT3bhcSJhJy1VC6k9tpTmvJMkhD20VfkRhJlBNy1AT1oKh5/fjSQl8OCnbwLgmTeOnzIVRamx\n0mA0yQydm9bpC/BytPm09VA9kVSrL+i+qy/i7V+X0dcbJjWYQk9KO5lt06h4+zqq918FgIS91OdU\nE0jws865jG1dR0lxZdBxtBlJS2J6VhJ1e8ph2QKe2VXJ0tnZw66vpNRYaJ/RJDN0fM/m0hqe21Nt\nHYxEKJphrWV9vKoLj8tJdprw8YVXEPD66U3pJMntIhxMIBxMoD25ja6MVhaEZtFwshNBSPYlEOru\nw5GXSWFbB4TCXHRzCevXLh5xfSWlxkKD0SQzNCCsKynkrtXW4pjzspN531VWjeXZ7WWkeJKobmtH\nmlJICiXSUFBJi6sDpydCc+5JGvOrSA2k4i9P5FhjA/NyCsjv7cfpcXHJ5XMIH6oiIy+Da9cuJS3J\no+tcq3OiwWiSGZgTtrm0hk5fgPQkD3estvqPjle38mqZ1XFdkOQiwZFEIOJny/4aMutn4O5PoK7w\nOEfn7KMzo4X0zhxWh5fiSugm0e2ht0VoON6IycvkQHkze18qZcaqRfz01WP6JE2dMw1Gk9DQZUSS\n0qzpG/QHWTR3Gi63k/bmbiIha35aS287iX0pFFYXM722iMzWacysLmZedxE7j1dQ39nBe5ZcTHpL\nFziEaQtnkFdVRzgQ4q6Hbh6siWmiRnUuNBhNQkOXEUlM8eJwCJfmpfKhVfPJykunt7UbtyOBzKQU\n/KaLiIngME6Se9PJbi3A25+EcfroCXVQUjibnXuaaK9sxjEjhxZfkIbNb8PMHOo8CYM1sU07K3Q+\nmjprGowmoYFlRAZqLA6Hg7TcdPrau+nyBZg2MxvT3WctspY1k/5QP639dQjWEABjDF3BFg43VeBx\nJHLTwkXktnTgTU7AUZjLjKYWaOvm6o+vGQxEj205AiLaga3OWlwf7YvIzcD3sNaxfswY880hxz8M\nfBFrbfhu4NPGmHeixyqj+8JAyBiz/AIW3fbSkzysLM7j3367l3l5qTjSEjl6pJ4NLx/icE8QjMF0\n9rDqmktZt6yAf/39n6jxleN1JhOKBAmZAEnONIpzZpHp66ehppWHHl5Lb1oyv/2bH0FmCotvunSw\njwrQ8UXqnMStZhSTovoWYDFwr4gMzXFTAVxnjLkE+DrR/GcxbjDGlGggGt6Glw+xt6KFX++oICk/\nC2+Pj3tXLeATd15GUnIC8xyGm5fNwkMK3/jgzbznkmIS3AaHOMn2TCcrIZ+LMpN49MevsHBhAXe+\n7zI6jpykt7yeObddyc2XzwY0W6w6P+JZMxpMUQ0gIgMpqg8NnGCMeSPm/B1Y+dHUGK1fu5hg+CDz\n8lI5fKyCk68fYufhOpJTvFx77UJeeXk/L7xexm/21VGYncw/3X0daXd4+Lff7mVvRQsmFGbH7/cQ\njBiuuHUZIlD169fxZiRz9YeuifftqUlmIqS3HvAJ4IWY1wZ4RUT2iMj6kS6azOmtR1OYk8I377uS\nT65dzHtuKYGIof1EA49tOULeklkkJnrY+/I+pmckUtPayw9fPEB6kocvvX8Zi/PTCB+qxt/Vx8Ib\nL+amK4rY+rtd7Hv1AAvvWc3Tu6rZtNPKk6ZP0dT5MCE6sEXkBqxg9MWY3auMMSVYzbyHROTa4a41\nxmwwxiw3xizPzc29AKW1p2WrFwLQfLCa+65dwB2rF/DQp2/i6NF6OFGPiRhy0xLp9AV46sX9lL+0\nF9PRg3NBAcf7Irx+oJYfff4JyE5DViyy3jSaikhXdVTnQzybaWNKUS0iS4HHgFuMMa0D+40xtdHf\nTSLyLFazb9u4lngCGC5vGkDuzGzSZ+aw68VSFn3gajaX1jBrdi6pc/OoPlADnnq21jTy6lOv42vt\nQjxunJfOpWB2DqsW5tP60h46qppZ882P8OH3LmXHsUZWFufxi61H8Qcj3Ld6vj5FU+fE7umtZwG/\nBe43xhyL2Z8MOIwx3dHtdcDXLljJbWzwMTt/maza6QuwaWcF+Ssvovt3Owj09rNxezlet4O+Gbk4\nE73Q0EZvdQuS4sUxMxfHjGwyMpJo6Ojj6L5q9n/n93DxHOatWkxhTgqFOSk888ZxNm4vB+CBNQu1\nA1udE7unt/4KkA38SKwmwcAj/Dzg2eg+F/BLY8yLcbgN24l9zA5WIPr6M7vZX92OycrEhMKkHj/J\nsqJs9la04nQIRYtncjw7DTGG4oIMMIZuf5CHbr6Y594oY+eXf4HH6+bub3yYlcV5g8uEaCoidT7Z\nPb31A8ADw1x3Arh03As4AQ08Zh+wubSG/dXtAExbOIPkZUX88QcvcOOPH2RvRSsRA0kJzsHgdOns\nLKpbeihr6KKyqQvvS7ugvo3P/8/nWfO+Ep554/gpNa/7r78oLvepJp8J0YGtzt66kkJyU605aDlp\nXj721bupK2/g4FPbWTQjA4BQ2DAvL40PXjWXo3Wd7CxvZllRNt2b97Dt6de55x8+QFtBLp2+gC4T\nosaNLq42yaUnebhuSQG/3nECl9PB4jVLmbtmKXuf2ML7rlpIanQO2+HaDpYVZbO/ug0TjtDzm9d4\n+rm3WPPh1bhvKOGxLUd4q6yJS+dkc9sVc7R/SJ13WjOaAu6+Zh4r5ueyv7qdzaU1/OPjD5Kal8Hv\nH/g+VwX7+PDq+dy3ej7z8tMxtS14fvkKR597i7u/cBt//8RDDKR73F/dxsZtZfoIX40LrRlNEiM9\n0oe/TJyNPf6j177O1+76Nt/92A9YcPlcMvPSaaprx5RW0u9xcf2X7+KTX7+bTl8AfzDMJbOymDMt\n9ZS5aEqdT1ozmiRON/BwuECVP2ca/7Llq1z1t+/F4XbS3tCBLxBGbrmS5T98kIe/dMfg+/56RwX7\nq9vIS0/k/uuKtYmmxoXWjCaJoY/0Y23aVcnGbWX4g+FTsoZsOVDHW9k5PPCtVdx19bxhg1bs4/vY\nx/oakNT5pjWjSeK0M+eNAcAfCJ0yh2xlcR4r5ueysjhvxPdIT/Jw//UXcf91xew41njKCpJKnU8a\njKaA21YU8cCahXg9rlOacjuONbKzvJkNLx8aU3AZuoKkUueTNtOmgIEaz0DA8QdC1LT04A+EWFaU\nw87yZr71XClfuL1k2JpVbPMttiNcqfNJg9EUkp7kwet28tiWIxyrtwY33rd6Pm6nDNZ2YvuOVhbn\nseNYI/7nVVp7AAAMVUlEQVRgmI3bygBr1PXAUrPad6TOJw1GU8xAjWZlcR5LZzcOvo6t7Qw8mdtX\n1ToYsGIzgHzruVJ2lltrQ2nmWHW+aDCaYE43nmik84BTrhkIIIU5KYPnxwaV4QLWwGc988ZxdpY3\ns2J+rjbV1HmlwWiCGW6JkNHOA0a8ZrjgNjCwcdPOisEF1Abo4vtqvGgwmmBON55otPOGuyY2aMX2\nBW0urRlcq8jrdg4GsaGrAih1vmgwmmDGGgxi01yvKyk85ZrY2tDK4jz2VbWysjjvXYFJ1ypSF5IG\no0lspCbd0CbczvLmUzqzB5pgulaRupA0GE1iIzXpRmrCaRNMxdNZj8AWka+cz4Ko82+kKSKxQUfH\nCym7OJfpIO9aDvZMicjNInJURMpF5JFhjouI/Ff0+D4RuWys16rRaYohZSenbaaJSNdIh4DEc/ng\nmPTWa7ESOO4SkU3GmEMxp90CLIj+XAk8Clw5xmvVKMb6ZE6pC2G0PqMO4ApjTOPQAyJyrv+cjpre\nOvr658YYA+wQkQwRmQ7MGcO1ahTaR6TsZLRm2s+B2SMc++U5fvZY0luPdM6YU2NP5fTWSk0kpw1G\nxph/MMbsFJGNIvJJEVkYc+yLp7vWLjS99ak6fYFT1jRSyi7G2oH9U2A68H0ROSEivxGRvzvHzx5L\neuuRzhlTamz1btpprexqTOOMjDGvisg24ArgBuBTwBLge+fw2aOmtwY2AQ9H+4SuBDqNMfUi0jyG\na9UwtNNa2dWYgpGIbAGSgTeB7Vid2k3n8sFjTG/9PHArUA74gI+f7tpzKc9UoZ3Wyq7GOgJ7H3A5\ncDHQCXSIyJvGmL5z+fAxpLc2wENjvVYpNXGNtZn2OQARSQU+BvwMyAcSxq1kSqkpZazNtIeB1Vi1\no0rgcazmmlJKnRdjbaZ5ge8Ae4wxoXEsj1JqihprM+1b410QpdTUpnnTlFK2oMFIKWULGoyUUrag\nwUgpZQsajJRStqDBSCllCxqMlFK2oMFIKWULGoyUUragwUgpZQsajJRStqDBSCllCxqMlFK2oMFI\nKWULcQlGIpIlIi+LSFn0d+Yw5xSKyKsickhEDsZmIxGRr4pIrYiURn9uvbB3oJQ63+JVM3oE2GKM\nWQBsib4eKgT8b2PMYmAl8JCILI45/l1jTEn0R9fCVmqCi1cwuh14Mrr9JHDH0BOMMfXGmLej293A\nYUbIGquUmvjiFYzyjDH10e0GIO90J4vIHGAZ8FbM7s+IyD4ReXy4Zl7MtZreWqkJYNyCkYi8IiIH\nhvm5Pfa8aDoic5r3SQF+A3zWGNMV3f0oMBcoAeqBb490vaa3VmpiGOuC/GfMGHPTSMdEpFFEpkez\nw04Hhk0IKSJurED0/4wxv41578aYc/4b+MP5K7lSKh7i1UzbBHw0uv1R4LmhJ4iIAD8FDhtjvjPk\n2PSYl3cCB8apnEqpCyReweibwFoRKQNuir5GRApEZODJ2DXA/cCNwzzC/w8R2S8i+4AbgM9d4PIr\npc6zcWumnY4xphVYM8z+OuDW6PZrgIxw/f3jWkCl1AWnI7CVUragwUgpZQsajJRStqDBSCllCxqM\nlFK2oMFIKWULGoyUUragwUgpZQsajJRStqDBSCllCxqMlFK2oMFIKWULGoyUUragwUgpZQsajJRS\ntqDBSCllCxqMlFK2oMFIKWULtk1vHT2vMrrWdamI7D7T65VSE4ed01sPuCGawnr5WV6vlJoAbJve\nepyvV0rZjN3TWxvgFRHZIyLrz+J6TW+t1AQxbqmKROQVIH+YQ1+OfWGMMSIyUnrrVcaYWhGZBrws\nIkeMMdvO4HqMMRuADQDLly8f8TylVHzZOr21MaY2+rtJRJ4FVgDbgDFdr5SaOOyc3jpZRFIHtoF1\n/CWN9ajXK6UmFjunt84DXhORd4CdwB+NMS+e7nql1MRl5/TWJ4BLz+R6pdTEpSOwlVK2oMFIKWUL\nGoyUUragwUgpZQsajJRStqDBSCllCxqMlFK2oMFIKWULGoyUUragwUgpZQsajJRStqDBSCllCxqM\nlFK2oMFIKWULGoyUUragwUgpZQsajJRStqDBSCllC7ZNby0iF0XTWg/8dInIZ6PHvioitTHHbr3w\nd6GUOp9sm97aGHM0mta6BLgc8AHPxpzy3YHjxpjnh16vlJpYJkp66zXAcWNM1biWSikVN3ZPbz3g\nHuCpIfs+IyL7ROTx4Zp5SqmJZdyCkYi8IiIHhvm5PfY8Y4wBRkw7LSIe4DbgmZjdjwJzgRKgHvj2\naa5fLyK7RWR3c3PzudySUmoc2Tq9ddQtwNvGmMaY9x7cFpH/Bv5wmnJsADYALF++fMSgp5SKL9um\nt45xL0OaaNEANuBO/pL2Wik1Qdk5vTUikgysBX475Pr/EJH9IrIPuAH43IUptlJqvNg2vXX0dS+Q\nPcx5949rAZVSF5yOwFZK2YIGI6WULWgwUkrZggYjpZQtaDBSStmCBiOllC1oMFJK2YIGI6WULWgw\nUkrZggYjpZQtaDBSStmCBiOllC1oMFJK2YIGI6WULWgwUkrZggYjpZQtaDBSStmCBiOllC1oMFJK\n2UJcgpGI3CUiB0UkIiLLT3PezSJyVETKReSRmP1ZIvKyiJRFf2sSR6UmuHjVjA4A7we2jXSCiDiB\nH2LlTVsM3Csii6OHHwG2GGMWAFuir5VSE1hcgpEx5rAx5ugop60Ayo0xJ4wxAeBpYCAb7e3Ak9Ht\nJ4E7xqekSqkLJS6pisZoBlAT8/okcGV0O88YUx/dbgDyRnoTEVkPrI++7BeRyZjwMQdoiXchxslk\nvbfJel8Xne2F4xaMROQVIH+YQ182xpwug+wZMcYYERkxbXVsemsR2W2MGbGPaqKarPcFk/feJvN9\nne214xaMjDE3neNb1AKFMa9nRvcBNIrIdGNMfTTVddM5fpZSKs7s/Gh/F7BARIpExAPcA2yKHtsE\nfDS6/VHgvNW0lFLxEa9H+3eKyEngKuCPIvJSdH+BiDwPYIwJAQ8DLwGHgV8ZYw5G3+KbwFoRKQNu\nir4eiw3n8TbsZLLeF0zee9P7GkKMGbG7RSmlLhg7N9OUUlOIBiOllC1M6mB0rtNO7Gqs02FEpFJE\n9otI6bk8ch1vo/39xfJf0eP7ROSyeJTzbIzh3q4Xkc7od1QqIl+JRznPlIg8LiJNI43bO6vvzBgz\naX+ARViDsLYCy0c4xwkcB+YCHuAdYHG8yz7Kff0H8Eh0+xHg30c4rxLIiXd5R7mXUf/+wK3AC4AA\nK4G34l3u83hv1wN/iHdZz+LergUuAw6McPyMv7NJXTMy5z7txK4m03SYsfz9bwd+biw7gIzo+DK7\nm4j/bY2JMWYb0HaaU874O5vUwWiMhpt2MiNOZRmrsU6HMcArIrInOi3Gjsby95+I3xGMvdxXR5sy\nL4jIkgtTtHF3xt+ZneemjcmFmnZyoZ3uvmJfGHPa6TCrjDG1IjINeFlEjkT/RVP28TYwyxjTIyK3\nAr8DFsS5THEx4YORGd9pJ3FzuvsSkTFNhzHG1EZ/N4nIs1jNBrsFo7H8/W35HY3BqOU2xnTFbD8v\nIj8SkRxjzESfRHvG35k2004/7cSuRp0OIyLJIpI6sA2sw1pHym7G8vffBHwk+oRmJdAZ00y1s1Hv\nTUTyRUSi2yuw/p9sveAlPf/O/DuLd6/8OPf434nVVu0HGoGXovsLgOeH9Pwfw3ry8eV4l3sM95WN\ntahcGfAKkDX0vrCe4LwT/Tlo5/sa7u8PfAr4VHRbsBbaOw7sZ4Qno3b8GcO9PRz9ft4BdgBXx7vM\nY7yvp4B6IBj9f+wT5/qd6XQQpZQtaDNNKWULGoyUUragwUgpZQsajJRStqDBSCllCxqMlO2JyOXR\n1QfKozPBJd5lUuefBiM1ETwKfBJrmsQC4Ob4FkeNBw1GKu5E5O9F5G+j298VkT9Ft28UkS1AmjFm\nh7EGxf2cib1KgRqBBiNlB9uB1dHt5UCKiLij+zZjjfAdMFFm7KszpMFI2cEe4HIRScOauvMmVlBa\nHd1WU8CEn7WvJj5jTFBEKoCPAW8A+4AbgPlY8+9mxpw+UWbsqzOkNSNlF9uBL2AtcbIda9LlXmPN\n9O4SkZXRp2gfQZN2TkoajJRdbAemA28aYxoBf3QfwIPAY0A51izwF+JSQjWudNa+UsoWtGaklLIF\nDUZKKVvQYKSUsgUNRkopW9BgpJSyBQ1GSilb0GCklLKF/w9PaE+/VUfCqQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113b41da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "w0, w1 = np.meshgrid(np.linspace(-1, 1, 100), np.linspace(-1, 1, 100))\n",
    "w_grid = np.array([w0, w1]).transpose(1, 2, 0)\n",
    "plt.contour(w0, w1, multivariate_normal.pdf(w_grid, mean=analytical_model.w_mean, cov=analytical_model.w_cov))\n",
    "w_sample = np.asarray(sample[\"w\"])\n",
    "plt.scatter(w_sample[:, 0], w_sample[:, 1], c=\"steelblue\", s=1)\n",
    "plt.xlabel(\"w0\")\n",
    "plt.ylabel(\"w1\")\n",
    "plt.xlim(-1, 1)\n",
    "plt.ylim(-1, 1)\n",
    "plt.gca().set_aspect(\"equal\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
